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ABSTRACT 

In addition to the Sun, six other stars are known to harbor multiple planets and debris disks: HD 
69830, HD 38529, HD 128311, HD 202206, HD 82943 and HR 8799. In this paper we set constraints 
on the location of the dust-producing planetesimals around the latter four systems. We use a radiative 
transfer model to analyze the spectral energy distributions of the dust disks (including two new Spitzer 
IRS spectra presented in this paper), and a dynamical model to assess the long-term stability of the 
planetesimals' orbits. As members of a small group of stars that show evidence of harboring a multiple 
planets and planetesimals, their study can help us learn about the diversity of planetary systems. 
Subject headings: circumstellar matter — Kuiper Belt — infrared: stars — planetary systems — stars: 
HD 128311, HD 202206, HR 8799, HD 82943 
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1. INTRODUCTION 

Surveys with the Spitzer Space Telescope have been 
spectacularly successful at identifying infrared excess 
emission associated with planetary debris disks around 
A- through K-type stars. The excesses at 70 ^m are 
associated with cool dust located at distances from the 
stars analogous to the position of the Kuiper Belt (KB) 
in the solar system (Moro-Martin et al. 2008), although 
much larger numbers of objects must lie in these exo- 
Kuiper belts to account for the detected level of emis- 
sion (Trilling et al. 2008; Hillenbrand et al. 2008; Car- 
penter et al. 2009). The most sensitive studies find 70 
/zm excesses around '^ 15-20% of mature solar- type stars 
over the entire 10 Myr to 10 Gyr age range (Trilling et 
al. 2008). The excesses at 24/im are generally associated 
with warmer dust and disappear relatively quickly as the 
host star ages; about 30% of the stars of the Pleiades age 
{^ 120 Myr) show excess emission (Sierchio et al. 2010) 
whereas by the age of Praesepe (~ 600 Myr), the 24 /im 
excesses have almost completely disappeared (Meyer et 
al. 2008; Caspar et al. 2009). Because the expected life- 
times of the debris dust grains are much shorter than the 
ages of the stars, it is inferred that the dust originates 
from coUisional activity in reservoirs of planetesimals left 
over from the planet formation process (hence the term 
debris dust). To sustain the dust production, it is neces- 
sary that large planetesimals (1000 km-sized) or unseen 
planets stir the planetesimals so they continue to collide 
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with each other. Planets are also responsible for con- 
straining the planetesimals in some zones and clearing 
them from others, thus determining much of the struc- 
ture of the debris system. The way the patterns of debris 
disks activity decay with age is consistent with the ex- 
pectation that the inner zones of a planetary system have 
relatively short dynamical time scales, whereas dynami- 
cal activity unfolds slowly at the distance of the Kuiper 
Beh. 

A highlight of the recent surveys is the first detection 
of debris disks around stars with planets (Beichman et 
al. 2005; Moro-Martin et al. 2007a). These planets orbit 
within several AU of their parent star, whereas the cold 
dust emitting the observed far-IR radiation generally re- 
sides tens of AU away. Despite the separation between 
the dust and planets, it is still possible for the planet 
to shape the structure of the dust (and planetesimals) 
disk. In the multiple-planet system HD 38529, for exam- 
ple, secular resonances excited by planets at 0.13 and 3.7 
AU create regions at tens of AUs that are unstable for 
orbiting planetesimals (Moro-Martin et al. 2007b). 

Several additional systems have been identified as hav- 
ing both multiple planets (capable of exciting secular res- 
onances) and orbiting debris (indicating the presence of 
planetesimals; Bryden et al. 2009, Su et al. 2009). Pre- 
sumably there are many more stars with debris-disk ex- 
cesses that also harbor multiple but undiscovered plan- 
ets. Detailed studies of the known examples can reveal 
aspects of their behavior that help us understand the di- 
versity of planetary systems. In this paper we study four 
of these systems: HD 128311, HD 202206, HD 82943 
and HR 8799. In Section 2 we describe the planet and 
debris dust detections for each one of these systems, and 
present new Spitzer IRS detections for the debris disks 
around HD 202206 and HD 82943. In Section 3 we use 
a radiative transfer model to identify the range of pa- 
rameters (dust mass and dust location) that would fit 
the observed spectral energy distribution (SED). Due to 
the high fractional luminosity, the grain-grain coUisional 
time-scale is shorter than the Poynting-Robertson (P-R) 
time-scale for all these systems and therefore we expect 
the dust to trace the location of the dust-producing plan- 



etesimals. In Section 4, we use a dynamical model to as- 
sess the long-term orbital stability of the putative dust- 
producing planetesimals, taking into account the effect 
of secular resonances. Putting together the results from 
the SED and dynamical analysis, in Section 5 we discuss 
the potential location of the dust-producing planetesimal 
belts. 

2. DEBRIS DISKS DETECTED IN MULTIPLE-PLANET 
SYSTEMS 

Including the solar system, 19 planetary systems are 
currently known to have both orbiting debris and planets 
(Bryden et al. 2009), and while the majority are single- 
planet systems, seven of them are known to harbor multi- 
ple planets. These systems are described in Figure 1 and 
Table 3. Figure 1 also shows the locations of the plan- 
etesimal belts derived either from work in the literature 
or in this paper. Our estimates arise from SED fitting 
together with dynamical models that study the effects 
imposed by the planets on the stability of the planetes- 
imals. In this section we describe the infrared excesses, 
the planet detections, and the dynamical models for the 
planet orbits. 

For these multi-planet systems, the planet-dust inter- 
action has previously been studied in the case of the solar 
system, HD 38529 (Moro-Martin et al. 2007b) and HD 
69830 (Lovis et al. 2006, Lisse et al. 2007). In this paper 
we study the remaining four systems: HD 128311, HD 
202206, HD 82943 and HR 8799. 

2.1. Observations of HD 128311 

HD 128311 is a KO star located at 16.57 pc, with Te// 
== 4965 K, M* = 0.84 Mq and a metaUicity of [Fe/H] 
== -0.04 (Saffe et al. 2008). King et al. (2003) identi- 
fied this star as a possible member of the UMa moving 
group, suggesting an age of ~ 500 Myr, consistent with 
the estimate of 390 to 410 Myr by Saffe et al. (2005). 
The chromospheric activity index reported by Gray et al. 
(2003) indicates an age of ~ 560 Myr, using the calibra- 
tion of Mamajek & Hillenbrand (2008). Barnes (2007) 
found an age of 350 Myr from gyrochronology. All of 
these determinations are consistent within their errors 
with an age of 500 Myr, which also lies at the lower end 
of the range from isochrone fitting (Valenti & Fischer 
2005). 

2.1.1. Debris dust detections 

HD 128311 was observed by a Spitzer Guaranteed 
Time Observation program specifically targeting planet- 
bearing stars. The 5piteer observations are shown in Fig- 
ure 2 and Table 1. The star has a strong excess at 70 
/im (Beichman et al. 2005; Trilling et al. 2008), and no 
excess was detected with MIPS at 24 /im (Trilling et al. 
2008). No excess was seen in the IRS spectrum at 5-35 
fim (Beichman ct al. 2006). [HD 128311 has a MIPS 
24 /im flux of 60±1.2 nijy (l-cr uncertainty), somewhat 
lower than but consistent with the flux expected from a 
K-band extrapolation (64 mJy)]. 

Excesses at 70 /xm are relatively common. Hillenbrand 
et al. (2008) and Carpenter et al. (2009) report them for 
6-10% of the 328 FGK stars in the FEPS sample, while 
the generally deeper (relative to the photospheric level) 
measurements of Trilling et al. (2008) find them in 16.4 
± 2.9 % of the FGK stars in their sample. However, HD 



128311 may be exceptional in not having an excess at 
33 /im, because only 3 out of 152 FGKM stars surveyed 
by Lawler et al. (2009) showed a 70 /im excess with no 
corresponding 33 /im. 

2.1.2. Planet detections 

Radial velocity monitoring of HD 128311 has led to 
the discovery of two planets. We found that the orbital 
solution in Vogt et al. (2005) is not stable. In this paper 
we have considered the following two orbital solutions 
(see planetary parameters in Table 2): 

• Fit Al is the solution in Butler et al. (2006), with 
the 2 planets near the 2:1 mean motion resonance 

(MMR). 

• Fit A2 (labeled as fit I in Gozdzlcwski & Konacki 
2006) corresponds to the best fit, co-planar solution 
lying at the border of an island of stability related 
to the corotation of apsides, with the planets also 
in the 2:1 MMR. 

Gozdzlcwski and Konacki (2006) suggested that the ra- 
dial velocity data can also be explained by two planets in 
a non-co-planar 1:1 MMR. In the present paper, we do 
not consider this somewhat exotic solution, as we assume 
that the protoplanetary flattened disk results in a system 
where the planets and the dust-producing planetesimals 
are on the same plane. 

Direct planet searches in the outer planetary system 
with the MMT (BiUer et al. 2007) and VLT/NACO 
(Eggenberger et al. 2007) have led to non-detections, 
however, they can set no constraints on the presence of 
planets smaller than 10 Mj„p. 

2.2. Observations of HD 202206 

HD 202206 is a metal-rich G6 V star located at 46.3 
pc, with Te// = 5765 K, M, = 1.15 M©, a metallicity 
of [Fe/H] ^ 0.37±0.07 and a stellar age of 5.6±1.2 Gyr 
(Udry et al. 2001); Saffe et al. 2005 estimates an age of 
4.2 Gyr. 

2.2.1. Debris dust detections 

HD 202206 was observed by a Spitzer Cycle 4 General 
Observer program targeting planet-bearing stars that 
were missed by the earlier Guaranteed Time Observa- 
tions. The Spitzer observations are shown in Figure 
2; Table 1 summarizes the (synthetic-)photometry. A 
strong excess was detected at 70 /im, a factor of 13 above 
the steUar photosphere (Bryden et al. 2009). In this pa- 
per we present new IRS observations which show emis- 
sion in excess of the stellar photosphere for A > 25 /im 
(sec also Dodson-Robinson 2010, in preparation). No 24 
/im observations were made. We use the new IRS spec- 
trum to calculate fluxes in four synthetic bands. 

2.2.2. Planet detections 

From radial velocity monitoring, Correia et al. (2005) 
concluded there are two planets in this system: a very 
massive 17.5 Mj„p inner body and a less massive 2.4 
Mjtjp outer planet. This system is particularly interest- 
ing because of the high mass of the innermost body (a 
brown dwarf). If it formed in the circumstellar proto- 
planetary disk, its existence would imply that such disks 
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Fig. 1. — Schematic representation of the seven planetary systems known to harbor multiple planets and dust-producing planetesinials. 
The stars are represented by the orange/yellow circles, with the stellar mass, luminosity and effective temperature labeled to the left. The 
size of the orange circle is proportional to the cube root of the stellar mass, while the size of the yellow circle is proportional to the stellar 
luminosity. The planets are represented by blue symbols with sizes proportional to the cube root of the planet mass. The thin blue lines 
extend from periastron to apoastron. For a given system, there is a range of planetary configurations that can fit the observations (see 
the planets' orbital elements in Table 2); many of these possible configurations would look very similar to each other under this schematic 
representation, with the exception of HD 82943 and HR 8799; for each of these stars we show two possible planetary configurations that 
significantly differ from each other. The inferred location of the dust-producing planetesimals are represented by the black and the green 
lines. Each line corresponds to a possible solution of a single-component disk, showing the degeneracy of the problem (exceptions are 
HR 8799, where warm and cold dusty disk components are inferred to exist from the SED, and the Sun, with the Asteroid Belt and the 
Kuiper Belt). The black lines correspond to solutions that assume that the dust is composed of 10 /jm-sized grains, while the green lines 
correspond to models that assume a grain size distribution (see Table 3 for details); the asteroid and Kuiper belts are shown in red, the 
only two planetesimal belts that have been directly detected. 



can be extremely massive, whereas if it formed like a 
stellar companion this system would be the only known 
example of a circumbinary planet orbit. In the latter sit- 
uation, Nelson (2003) predicted that the interaction of 
the outermost planet with the viscous circumbinary disk 
could have been responsible for its inward migration and 
resonance trapping. 

From stability considerations, Correia et al. (2005) 
suggested that the system is in an island of stability 
around the 5:1 mean motion resonance (MMR), with 
Msini of 17.428 Mj„p and 2.436 Mj„p, semi-major axes 
of 0.830 AU and 2.542 AU, and eccentricities of 0.435 and 



0.267, for HD 202206 b and c, respectively. This solution 
has a {xu^y^^ of 1.67 and is stable over a 5 Gyr timescale. 
However, using n-body analysis that takes into account 
stability considerations and is well suited for multi-planet 
systems in low-order MMRs, Gozdziewski et al. (2006) 
concluded that the solution by Correia et al. (2005) is 
a local minimum. Assuming a co-planar system, their 
best fit (with r.m.s = 9.98 m/s and (xu^Y^^ = 1-53) has 
Msini of 17.624 Mj^p and 2.421 Mj„p, semi-major axes 
of 0.831 AU and 2.701 AU, and eccentricities of 0.433 and 
0.255, for HD 202206 b and c, respectively. [Gozdziewski 
et al. (2006) noted that an unexpected non-co-planar 



TABLE 1 
Spitzer Photometry 



Star 



13.20 ' 



19 



23.68 



Wavelength (/.tin) 
25 " 



32.5 



71.42'' 



155.9 



HD 128311 
HD 202206 
HD 82943 


192.45 ± 15.62 

57.57±3.80 
201.10±10.87 


28.94±1.95 
98.29±6.99 


60.0 ± 1.2 

N/A 
66.0±1.3 


17.22±1.16 


32.95 ± 4.82 
14.23±0.83 " 
46.98±3.27 


23.5 ± 3.3 " 

32.0 ± 3.1 " 

133±5'= 


50.0 ± 35.0 

N/A 
N/A 




Note. — Fluxes are given in mjy. The uncertainties are 1-a. 

^Synthetic fluxes calculated from the IRS spectrum using a square filter over the following wavelength ranges: 
13.2 fim (12.4-14 fim), 19 fim (18-20 /tm), 25 /tm (24-26 /tm), and 32.5 fim (30-35 ^im). The uncertainties in this 
case are calculated from the standard deviation of the fiux within the wavelength bin. 

''The 70 /tin fiux is color corrected assuming a dust temperature of 50 K and a correction factor of 0.893 (from the 
Spitzer Users Manual) . 

'^Measured fiux is significantly above that expected from the stellar photosphere. 

solution is close to the center of libration and has 
{x,/y/^ = 1.553. 

The radial velocity data acquired so far cannot exclude 
the presence of a Neptune-sized planet between 0.06 AU- 
0.12 AU, or a 10 M® planet between 0.02 AU-0.12 AU, 
or a 0.5 Mj„p mass planet outside 6.5 AU. No planetary 
companions have been found by direct imaging (a source 
detected by Chauvin et al. 2006 was a background ob- 
ject). 

2.3. Observations of HD 82943 

HD 82943 is a GO V star located at 27.46 pc, with Tg// 
= 5989 K, M* = 1.15 M© and [Fe/H] = 0.26 (Sousa et al. 
2008). Consistent measurements of the chromospheric 
activity have been reported by Wright et al. (2004) , Gray 
et al. (2006), and Saffe et al. (2005). Using the calibra- 
tion of Mamajek and Hillenbrand (2008), they imply an 
age of 5 Gyr. 

2.3.1. Debris dust detections 

HD 82943 was observed by a Spitzer Guaranteed Time 
Observation program targeting nearby late F, G, and 
early K stars. The Spitzer observations are shown in Fig- 
ure 2; Table 1 summarizes the (synthetic-)photonietry. 
Trilling et al. (2008) found a strong excess emission at 
70 /im and no excess at 24 /xm, deriving a minimum char- 
acteristic dust temperature of 69 K (corresponding to 22 
AU if assuming blackbody grains), and a fractional lu- 
minosity of L^Mst/Lstar ~ 10""*. In this paper we present 
new IRS observations that show how the spectrum rises 
above the stellar photosphere for A > 26 /im. 

2.3.2. Planet detections 

From radial velocity observations. Mayor et al. (2004) 
announced the presence of two planets around HD 82943. 
Ferraz-Mello et al. (2005) found that the orbital solution 
in the discovery paper is unstable and that the radial ve- 
locity observations can be fitted with a stable, co-planar 
solution in which the planets are locked in the 2:1 MMR. 
Gozdziewski & Konacki (2006) revisited the orbital solu- 
tion modehng data from CORALIE (Mayor et al. 2004) 
and from Keck-HIRES (Lee et al. 2006), confirming that 
the co-planar 2:1 MMR solution fits the observations. 
They also found a stable, non-co-planar solution in a 1:1 
MMR, but as we mentioned above in this paper we only 
consider co-planar configurations. The dynamical maps 



Fig. 2.— Spitzer observations of HD 128311, HD 202206 and 
HD 82943. The IRS spectrum is shown as a continuous line, while 
squares are the MIPS photometric points; in both cases the error 
bars correspond to 1-cr uncertainties. The dotted line is the stellar 
photosphere (approximated as a blackbody). 

solution is also possible, where the eccentricity of the in- 
nermost largest planet varies with larger amplitude than 
that of the outermost planet, while the inclination of the 
latter can assume almost any value]. In a more recent 
paper, Couetdic et al. (2009) used frequency map anal- 
ysis and updated radial velocity data to explore the long 
term stability of a wide range of possible orbital solu- 
tions. Adopting an updated stellar mass of 1.044 M©, 
they found acceptable co-planar configurations for incli- 
nations between 30°-90° with respect to the line of sight 
and favored an edge-on co-planar solution in which the 
bodies are in the 5:1 MMR (with (x,,^)^^^ = 1.4136). 
This solution has Msinj of 16.59 Mj„p and 2.179 Mj„p, 
semi-major axes of 0.8053 AU and 2.49 AU, and eccen- 
tricities of 0.431 and 0.104, for HD 202206 b and c, re- 
spectively. It is stable for more than 5 Gyr and differs 
from the solutions above in that the eccentricity of planet 
c is lower, allowing more stable regions outside the res- 
onances. Couetdic et al. (2009) explored the stability 
of test particles in the above configuration identifying 
two possible niches for stability: a < 0.12 and a > 6.5 
AU (for an integration time of 16000 years). Couetdic 
et al. (2009) noted that the solution above has a high 
resonant mode amplitude that will likely be dampened 
by dissipative forces. 

• Fit Bl in Table 2 is the edge-on co-planar con- 
figuration favored by Couetdic et al. (2009); this 



in Gozdziewski & Konacki (2006) revealed that near their 
unstable, best fit solution there are two narrow islands 
of stability associated with the 2:1 MMR (see orbital el- 
ements in Table 2): 

• Fit CI (labeled fit V in Gozdziewski & Konacki 
2006) lies in one of the two islands and corresponds 
to their best fit, rigorously stable, two-planet so- 
lution, characterized by co-rotation of the apsidal 
lines. [We found that the second solution they iden- 
tify in their Figure 9 becomes unstable after 27.8 
Myr and therefore we do not consider it in this pa- 
per; this solution lies at the very edge of the second 
island of stability of the 2:1 MMR, where slight dif- 
ferences in the initial conditions can result in a very 
different dynamical evolution.] 

Gozdziewski & Konacki (2006) and Beauge et al. (2008) 
discussed the inadequacy of the 2-body solution and 
speculated on the presence of a third planet to account 
for the large r.m.s. and the lack of convergence of some of 
the orbital elements (eccentricity and longitude of periap- 
sis) as the number of radial velocity data points increases. 
In this paper, we have also considered the following sta- 
ble, 3-planet, co-planar solutions (see orbital elements in 
Table 2): 

• Fit C2 is the solution in Gozdziewski & Konacki 
(2006), in which planets b and c are in the 2:1 MMR 
while the outermost planet d is in a low eccentricity 
non- resonant orbit. 

• Fit C3 is the solution in Beauge et al. (2008), in 
which the planets are in a Laplace ld:2c:4b MMR 
in a double asymmetric apsidal corotation reso- 
nance. 

2.4. Observations of HR 8799 

HR 8799 is a metal-poor A5 V star located at 39.4 pc, 
with Te// = 7500 K, M* = 1.5 M© and a metaUicity of 
[Fe/H] = -0.55. Marois et al. (2008) estimated a steUar 
age of 30-160 Myr; Moya et al. (2010) argued that the 
age is still unconstrained, with astroseismology analysis 
favoring an older age of ^ 1 Gyr, but requiring a better 
determination of the rotation velocity of the star. 

2.4.1. Debris dust detections 

HR 8799 was observed by Spitzer under programs 
50175 and 530; images and an SED analysis were pre- 
sented in Su et al. (2009). Deep IRS and MIPS obser- 
vations (including the MIPS SED mode from 55-95 /im 
- a wavelength range critical to constrain the shape of 
the SED), revealed a spatially resolved disk at 24 and 70 
/xm, whose outer boundary can be traced beyond 1000 
AU. The excess dust emission can be detected in the 
IRS spectra starting bellow 20 /im (less deep IRS obser- 
vations are discussed in Chen et al. 2009) and also at 
160 iim. Unresolved debris disk detections exist at 60 
and 90 /im (with ISO; Moor et al. 2006), at 850 /im 
(with JCMT/SCUBA; WiUiams & Andrews 2006) and 
at 1.2 mm (Sylvester et al. 1996). Preliminary analy- 
sis shows that the debris disk is spatially resolved with 
APEX/LABOCA at 870 /xm (Kalas et al. 2010). 



2.4.2. Planet detections 

Three distant planets have been detected by direct 
imaging at projected separations of 68 AU, 38 AU and 
24 AU (for HR 8799 b, c and d, respectively - Marois et 
al. 2008; Lafreniere et al. 2009); the planets' long orbital 
periods and the short baseline of the observations make 
the determination of the orbital parameters still uncer- 
tain. It is possible to estimate the masses of the directly 
detected planets from their luminosities by using evolu- 
tionary cooling models (assuming the age of the system 
is known). However, current models are discrepant and, 
in the case of HR 8799, the age of the system is un- 
certain, making the planets' masses poorly constrained; 
lower limits to the masses of HR 8799 b, c and d are in 
the ranges of 5-11, 7-13 and 7-13 Mj„p, with nominal 
masses of 7, 10 and 10 Mj„p, respectively (Marois et al. 
2008; Reidemeister et al. 2009). 

Fabrycky & Murray-Clay (2010) found that many or- 
bital solutions are unstable, including the face-on cir- 
cular configuration (for any range of reasonable planet 
masses). They tried to determine the system parame- 
ters by fitting simultaneously the astrometry (positions 
and proper motions), the planets' luminosities (that con- 
strain the range of possible planets' masses), and the re- 
quirement that the system is stable for at least 30-160 
Myr. 

• Fit Dl is the face-on, co-planar solution identified 
by Fabrycky & Murray-Clay (2010) that is stable 
for > 160 Myr even if the planets' masses are up to 
1.9 times their nominal values; in this solution the 
planets are locked into a Laplace ld:2c:4b MMR, 
allowing them to avoid close encounters (see orbital 
elements in Table 2). 

Reidemeister et al. (2009) have also explored the sta- 
bility of the HR 8799 system in the case of co-planar 
and initially circular orbits, where the orbital plane is 
allowed to have a range of inclinations (/ = 0-40°) and 
orientations (i7' — 0-180°). Here, / is the angle between 
the angular momentum vector and the vector toward the 
observer, and H.' is the angle between the North direction 
and the line of nodes (measured toward the East - see 
Figure 1 in Reidemeister et al. 2009). They integrated 
a grid of models over 100 Myr and identified as stable 
the configurations where the planets did not suffer close 
encounters during that period; these configurations (all 
co-planar) have / > 20° and $1' = 0-50°. We take into 
account the following considerations: 

• Figure 8 in Reidemeister et al. (2009), showing the 
models that fulfill the stability criteria (i.e. no close 
encounters during 100 Myr). 

• Figure 9 in Reidemeister et al. (2009), showing the 
models that fit the observed astrometry. 

• Inclination constrains of the equatorial plane of the 
star derived from the comparison of the measured 
rotational velocity w-sim and the expected rota- 
tional velocity for an A5 star: / = 7-22° (Lafre- 
niere et al. 2009). 

• Inclination constraints from astrometry: (1) as- 
suming circular, co-planar orbits and fixing the 



stellar mass to 1.5 M©, but letting the orientation 
of the orbit float: / -- 20° (Fabrycky & Murray- 
Clay 2010); (2) from observations with a 10 year 
baseline (using NICMOS archival data where the 
planets are detected): / = 13-23° (Lafreniere et 
al. 2009). 

• Inclination constraints from the debris disk mor- 
phology: (1) from Spitzer spatially resolved obser- 
vations: / < 25° (Su et al. 2009); (2) from pre- 
liminary analysis of 870 /im observations: / > 20° 
(Paul Kalas, private comm.). 

Based on the above considerations, we have adopted the 
following two additional non-face-on, co-planar configu- 
rations (see orbital elements in Table 2): 

• Fit D2 with (/ = 25°, ri' = 20°). 

• Fit D3 with (/ = 20°, n' = 45°). 

As for the planets' masses, we initially considered: (1) 
the nominal masses of 7, 10 and 10 Mj„p for HR 8799 
b, c and d, respectively; and (2) the low mass case with 
5, 7 and 7 Mj„p for HR 8799 b, c and d, respectively. 
However, we found that when adopting nominal planet 
masses, the configurations D2 and D3 were unstable. 

Based on the location of HR 8799 b in a pre-discovery 
2007 image (Metchev et al. 2009), Fabrycky & Murray- 
Clay (2010) ruled out stable co-planar configurations in 
which the innermost planet is very eccentric (e > 0.95) 
or retrograde. They also found very non-co-planar (non- 
resonant) configurations that are stable but we will not 
consider those because: (1) they imply we are observing 
the system at a special time; and (2) it is more reasonable 
to assume that the planets and the planetesimals that 
produce the dust are co-planar because they formed out 
of a flattened disk. 

In a recent study of HR 8799, Gozdziewski & Mi- 
gaszewski (2009) used a self-consistent, n-body analysis 
that takes into account stability considerations to carry 
out a quasi-global search for stable, co-planar planetary 
configurations that can fit the astrometric data and the 
astrophysical mass constraints for the star and the plan- 
ets. In their study, the planets and stellar masses are 
free parameters and are allowed to vary within their 1- 
a range, with M^ = 7^2 Mj„p, M^ = 10±3 Mj„p, M^ 
= 10±3 Mjup and M* = 1.5±0.3 Mq. As in Fabrycky 
& Murray-Clay (2010) and Reidemeister et al. (2009), 
the goal of their study is to use the requirement of sta- 
bility to identify long-lived configurations that may lie 
close to unstable best-fitting solutions, mitigating the er- 
rors introduce by the short baseline of the observations. 
Gozdziewski & Migaszewski (2009) found two long-term, 
stable, non-face-on, co-planar solutions (see orbital ele- 
ments in Table 2): 

• Fit D4 (/ == 15.5°, O' = 11.2°): inside a smah 
island of regular motion; the planets are in the 
ld:2c:4b MMR; the solution becomes unstable af- 
ter 400 Myr. 

• Fit D5 (/ == 11.4°, n' = 357.2°): also inside a 
small island of stable motion; the inner two planets 
are in the lc:ld MMR, and the outermost planet 
is in a narrow island close to the (lc:ld):3b MMR; 



the inner planets are stable for > 3 Gyr, while the 
3-planet system has a regular motion only during 
the first 600 Myr. 

Table 2 lists the orbital elements of the five planetary 
configurations (Fits Dl -D5) we have adopted for our test 
particle simulations of HR 8799 in Section 4.4. Our nu- 
merical runs last 160 Myr (the upper limit to stellar age 
estimated by Marois et al. 2008); there is no guarantee 
that the systems will remain stable on longer time-scales. 

3. SED MODELING 

We modeled the excess dust emission using the ra- 
diative transfer code developed by Wolf & Hillenbrand 
(2003). For a study on how the SEDs depend on the 
grain composition we refer to Wolf & Hillenbrand (2003) 
and Moro-Martin, Wolf & Malhotra (2005). For the grain 
properties, we assumed a silicate composition with opti- 
cal constants from Weingartner & Draine (2001), and 
two assumptions for the particle radius (&): 

• A single grain size of b = 10 fim. This size was 
chosen because these grains emit efficiently at 70 
/im. We favor this scenario because it represents 
large grains in general, where the grains would be 
located at the blackbody equilibrium distance from 
the central star (modulo bulk optical properties). 
Given the lack of spectral features in most debris 
disks (see e.g. Beichman et al. 2006 and Lawler et 
al. 2009), it is inferred that small grains with super- 
thermal behavior are mostly absent, therefore, this 
single grain size scenario provides a plausible lim- 
iting case for the placement of the grains. 

• A particle size distribution following a power law. 
The placements of the emitting zones are rather 
uncertain. In addition to the single grain size 
scenario above, the other bounding case would 
be given by assuming a power law down to the 
blow-out size. We used n{b) oc b~^'^ , with bmax 
= 10 ^im and b,nin — buow The blow-out size, 
bbiow, is the grain radius for which /3, the ratio 
between the radiation pressure force and the gravi- 
tational force, equals 1/2; for spherical grains, /3 — 
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Lamy & Soter, 1979), where p and b are the den- 
sity and radius of the grain in cgs units and Qp^ is 
the radiation pressure coefficient, a measure of the 
fractional amount of energy scattered and/or ab- 
sorbed by the grain. Qpr is a function of the phys- 
ical properties of the grain and the wavelength of 
the incoming radiation; the value we use, < Qpr >, 
is an average integrated over the stellar spectrum. 

The presence of larger and colder grains is not con- 
strained by the Spitzer observations, so generally these 
models lead to lower limits of the dust disk mass. 

We now discuss our assumptions for the disk geometry. 
We adopted a constant surface density disk (E oc r"). 
The dust disk outer radius. Rout, cannot be determined 
from the Spitzer data alone. Given that the debris disks 
around HD 128311, HD 202206 and HD 82943 are not 
spatially resolved at 70 ^im (with a PSF FWHM of 18"), 
and that the debris disk around e Eri (of similar spectral 
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Planetary and stellar parameters for the dynamical simulations 
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Note. — Planetary parameters used in the dynamical simulations in Section 4. a and e are the semi-major 
axis and eccentricity of the planet; in all cases the orbits considered are co-planar (i = 0); w is the longitude 
of periastron; Q is the longitude of the ascending node; M is the mean anomaly. The orientation of the orbit 
is given by angles I and Q': I is the inclination of the orbital plane with respect to the plane of the sky, and 
r2' is the angle between the North direction and the line of nodes (measured toward the East — sec Fig. 1 in 
Reidemeistcr et al. 2009). References are: B06: Butler et al. (2006); G06: Gozdziewski & Konacki (2006); C09: 
Couetdic et al. (2009); FIO: Fabrycky & Murray-Clay (2010); B08: Beauge et al. (2008); R09: Reidemeistcr et 
al. (2009); G09 Gozdziewski & Migaszewski (2009). 
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Fig. 3.— Observed and modeled SEDs for HD 128311. The dot- 
ted line is the stellar photosphere. The Spitzer photometric points 
{MIPS and synthetic photometry from IRS) are represented by 
squares with l-cr error bars. The continuous lines correspond to 
the modeled SEDs and include the emission from the photosphere 
and from a disk of dust composed of particles with optical proper- 
ties typical of astronomical silicates. Panels (a), (b) and (c) corre- 
spond to 10 fim size grains. Panel (d) corresponds to a grain size 
distribution given by n(b) oc b~^'^ and a maximum grain radius of 
bmax = 10 /xm; regarding the minimum grain radius, we adopt the 
arbitrary value of femin = 0.17 /^m because the low luminosity of 
HD 128311 results in the absence of a blow-out size. The dust disk 
has a mass yidust ^nd extends from Ri„ to Rout with a constant 
surface density. Ri„ and yidust a^r^ the free parameters; Rout is 
kept fixed at 10 AU, 50 AU and 100 AU (indicated at the top of 
each panel). The models in yellow are those with a x^ probability 
P(X^ I >^) < 0.683; light blue for P(x^ | v) > 0.683, i.e. models 
that can be excluded with l-tr certainty; and dark blue for P(x^ 
I v) > 0.9973, i.e. models that are excluded with 3-cr certainty. 
The triangles correspond to the 5-(T/lhr sensitivity limits for Her- 
schel/PACS: 3.75 mjy (60 fim-85 fim), 4.1 mjy (85 fim-130 fim) 
and 5.75 mjy (130 ^m-210 ^m). 

type) was well constrained by the Spitzer data to be ~100 
AU (Backman et al. 2009), we considered three possible 
outer disk radii: Rout = 10 AU, 50 AU and 100 AU. To 
explore different belt widths, the inner disk radius. Km, 
is allowed to vary from the sublimation radius (where 
T,„b = 1550 K) to Rout. 

We further assumed that the dust disk is optically thin 
(which is supported by the low fractional luminosities ob- 
served for HD 128311, HD 202206 and HD 82943, with 
^dust/^star < 10^"* ~ see Table 3). Under this scenario, 
only scattering, absorption and reemission of stellar radi- 
ation by dust grains were taken into account (neglecting 
scattering and dust heating from the dust infrared radi- 
ation). With the parameters described above, we calcu- 
lated the dust disk emission at 200 wavelengths from 3 
fiia to 600 /im. 

3.1. HD 128311 

Beichman et al. (2006) modeled the IRS and MIPS 70 
/im data with a single population of 10 /zm amorphous 
sihcate grains located at Rdust > 15 AU [Tdust < 50 K), 
while the modeling of the MIPS 24 /im and 70 /im data 
in Trilhng et al. (2008) implied Rdust > 5.1 AU (T^^^t < 
106 K), with a fractional luminosity ^dust/^star — 1-3- 
2.7-10-^ 

To model the dust emission, the stellar contribution 
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Fig. 4. — Parameter space of the modeled SEDs in Figure 3. Each 
point of these two-dimensional grids represents a modeled SED 
from Figure 3, where Ri„ and M^dust a-re the two free parameters. 
The models in yellow are those with a x^ probability P(x^ I '^) 
< 0.683; light blue for P(x^ I J^) > 0.683, i.e. models that can be 
excluded with 1-cr certainty; and dark blue for P(x'^ I J^) > 0.9973, 
i.e. models that are excluded with 3-cr certainty. 

needs to be subtracted from the observed SED. Because 
the slope of the ZRS' spectra of HD 128311 indicates there 
is no dust excess at A< 34 /xm, we pin the photosphere to 
the MIPS 24 /j,m flux; for simplicity, the stellar emission 
was modeled as a blackbody with Tg// — 4965 K, L* = 
0.24 L0 and a distance of 16.57 pc. 

For a stellar mass of M* = 0.84 Mq, a grain density 
p = 2.5 g cm~^, and adopting the optical constants of 
astronomical silicates, we find that the maximum /3- value 
is 0.26, i.e. there is no blow-out size because /3 < 0.5. 

Figure 3 shows the SED models compared to the SED 
observed by Spitzer. We calculated the synthetic pho- 
tometry of each SED model using the MIPS and IRS 
filter profiles from the Spitzer Users Manual, with effec- 
tive wavelengths at 13.2 /im, 23.68 /im, 32.5 //m, 71.42 
/im and 155.9 /im. Then, for each model we calculated 
X^ and the x^ probability, where the fit is done to the 
5 photometric points in Table 1 with two free parame- 
ters (Rmin and M^dust)- [Becausc the uncertainty in the 
160 /im flux is large, even though we are including it 
in the SED fit, in practice its contribution is negligible]. 
The colors indicate the goodness of the fit. Figure 4 
shows the parameter space of the SED models in Figure 
3. The yellow area corresponds to combinations of Ri„ 
and Mdust that lead to a x^ probability of P(x^ I ^) < 
0.683, light blue is for 0.683 < P(x^ | i^) < 0.9973, and 
dark blue if for P(x^ | i^) > 0.9973. The modeling of 
the SED is degenerate. For the single grain size models, 
valid SED fits can be obtained for: Rout = 10 AU and 
R„ > 5 AU; Rout = 50 AU and R^ > 2 AU; and Rout 
= 100 AU and Ri„ > 2 AU. In all three cases, narrow 
belts with widths of 10% the disk radius could also fit 
the observed SEDs. Here, we would like to emphasize 
the lower limits to Ri„ because of two reasons: (1) For 
HD 128311, HD 202206 and HD 82943 the dust is located 
outside the orbits of the planets and, therefore, it is at 
the disk inner edge where the gravitational effects of the 
planets are stronger (i.e. it is at that location that there 
is a closer correspondence between the SED models and 
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Fig. 5.— Same as Figures 3 and 4 but for HD 202206. Of the 
range of models explored, only the two sets of models shown in the 
panels could fit the observations. 

the dynamical models). (2) The narrowest belts might 
only be justified in the presence of additional planetary 
perturbers for which we have no evidence so far. 

As we mentioned before, we favor the single grain size 
models above because they represent large grains in gen- 
eral, where the grains would be located at the blackbody 
equilibrium distance from the central star. However, for 
completeness and to explore the other limiting case we 
have also considered the case of a grain size distribution. 
For these models, valid fits can be obtained for Rom* = 
100 AU and Ri„ > 52 AU. A small disk is ruled out for 
this grain size distribution. We can exclude the presence 
of a significant population of small grains unless the disk 
is large, 100 AU, in which case there would be a depletion 
of small grains inside 52 AU. 

3.2. HD 202206 

The slope of the IRS spectrum of HD 202206 indicates 
there is no dust excess at A< 25 /im; therefore, we used 
the IRS synthetic photometry at 15-19 /Ltm to determine 
the photospheric emission level. The stellar emission was 
modeled as a blackbody with Teff — 5764 K, L, = 0.92 
Lq and a distance of 46.3 pc. 

For a stellar mass of M* = 1.044 M0, a grain density 
p = 2.5 g cm~^ and adopting the optical constants of 
astronomical silicates, the blow-out size is hbiow = 0.5 
/im. 

For the SED modeling, we followed the same scheme 
as described above for HD 128311. Figure 5 (top) shows 
some of the SED models computed; over-plotted are the 
Spitzer observations. For each model, we calculated the 
synthetic photometry using the MIPS and adopted IRS 
filter profiles at effective wavelengths 13.2 /im, 19 /im, 25 
/im, 32.5 /im and 71.42 /im. We then calculated x^ and 
the x^ probability distribution, where the fit is done to 
the observed 5 photometric points in Table 1 with two 
free parameters (Rmin and Mdust)- The colors indicate 
the goodness of the fit. Of the wide range of models ex- 
plored, only the two sets of models shown can fit the ob- 
servations. Figure 5 (bottom) shows the parameter space 
(Rmm and yidust) of the models in the top two panels. 
The SED modeling is degenerate, but compared to HD 
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Fig. 6.— Same as Figures 3 and 4 but for HD 82943. Of the 
range of models explored, only the two sets of models shown in the 
panels could fit the observations. 

128311, the case of HD 202206 is better constrained be- 
cause of the detection of an excess starting near 25 /^m. 
For the models that assume a single grain size, the best 
fits are given by: Rout = 50 AU and 10 AU < Ri„ < 
20 AU; and Rout = 100 AU and 6AU < Ri„ < 10 AU. 
For this case we can rule out the presence of a small 10 
AU disk. We can also rule out the models that assume 
a distribution of grain sizes with a power-law index of 
-3.5, brmn = hiow and bmax = 10 /im. The latter might 
not surprising because collisional disk models by several 
authors (e.g. Krivov et al. 2006; Thebault & Augereau 
2007; Mueller et al. 2010, Krivov 2010) conclude that 
the size distribution shows substantial deviations from a 
power-law near the blowout size regime, suggesting that 
bmin would need to be two or three times larger than 
bbiow, and that the power-law index may be smaller than 
3.5, to mimic a secondary maximum seen in the size dis- 
tribution at hundreds of microns. 

3.3. HD 82943 

Because the IRS spectrum of HD 82943 indicates there 
is no dust excess at A< 22 /im, we used the IRS syn- 
thetic photometry flux calculated between 15-19 /im to 
pin down the photosphere. For simplicity, the stellar 
emission was modeled as a blackbody with Tg// — 5989 
K, L* — 1.25 Lq, and a distance of 27.46 pc. 

For a stellar mass of M* = 1.15 Mq, a grain density 
p = 2.5 g cm~'^ and adopting the optical constants of 
astronomical silicates, the blow-out size is bbiow = 0.6 
/im. 

Figure 6 (top) shows the SED models over-plotted on 
the Spitzer observations. For each model, we calculated 
the synthetic photometry using the MIPS and adopted 
IRS filter profiles at effective wavelengths of 13.2 /im, 19 
/im, 23.68 /im, 32.5 /im and 71.42 /im. We then calcu- 
lated x^ and the x^ probability distribution, where the 
fit is done to the observed 5 photometric points in Table 
1 with two free parameters (Rmm, ^dust)- The colors in- 
dicate the goodness of the fit. Figure 6 (bottom) shows 
the parameter space of the models in the top panels. 

For the models that assume a single grain radius of 10 
/im, the best SED fits are given by: Kout = 50 AU and 
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16 AU < R,„ < 44 AU; and Rout = 100 AU and 12 AU < 
Ri„ < 26 AU. For this case we can rule out the presence 
of a small 10 AU disk. As for HD 202206, we can also rule 
out models that assume a grain size distribution with a 
power-law index of -3.5, bmin — bhiow and bmax ~ 10 /^m. 

3.4. HR 8799 

The Spitzer observations of HR 8799 and the analysis 
of its SED and spatially resolved images were presented 
in Su et al. (2009). Here we summarize the main results. 
The stellar emission is best fitted by a Kurucz model with 
Te// = 7500 K, R, = 1.4 Rq, log(g) = 4.5 and sub-solar 
abundances, with a resulting L* = 5.7 Lq (consistent 
with a young age). The stellar mass is M* = 1.5 Mq. 
As in Sections 3.1 and 3.2, Su et al. (2009) assumed 
spherical grains with optical properties characteristic of 
astronomical silicates, a bulk density of 2.5 g cm~'^, and 
a size distribution following n{b) oc b~^'^; the blow-out 
size for HR 8799 is 2 fj,m. Su et al. (2009) found that 
the stellar-subtracted images and SED (after applying 
color-corrections) can be best fit by the following multi- 
component disk model: 

• Unresolved warm disk with a characteristic temper- 
ature of r^ 150 K. 
This component was modeled with a flat surface 



density, E 



The maximum inner radius is YL, 



6 AU (based on temperature arguments) and the 
outer radius is Rout = 15 AU. The grain sizes range 
from 1 to 4.5 /xm (with the spectral shape between 
10 and 20 /im requiring the presence of /zm-sized 
grains and the absence of sub-micron grains). If we 
were to assume all the grains are bound (> 2 /im), 
the outer radius of this warm component would be 
Rout ^10 AU (instead of 15 AU). The estimated 
dust mass is 1.1-10"^ M0. This is a lower limit be- 
cause it is not possible to constrain the mass con- 
tributed by larger grains as there are no long wave- 
length detections of this warm component. This 
component dominates the emission at IRS wave- 
lengths and the unresolved emission at 24 /im (it is 
too warm to contribute to the unresolved emission 
at 70 /im). The upper limit to the dust mass < 6 
AU is 2.3-10~^ Mg (assuming 10 /im grains). 

• Unresolved cold disk with a characteristic tempera- 
ture of '^45 K. 

This component was also modeled with a flat sur- 
face density, E '--^ r°. The inner radius, Ri„ = 90 
AU, is not constrained from imaging but from the 
SED: derived from the characteristic temperature 
of the cold component assuming blackbody grains 
10 /im in size. The sharp inflection between the IRS 
and the MZPS'-SED spectra indicates that there is 
a sharp inner edge of the cold dust component (i.e. 
there is very little dust located between the warm 
and the cold dust components). The outer radius is 
set to RoMf = 300 AU, but is less well constrained 
because of the lack of spatially resolved sub-mm 
observations; the value of 300 AU is chosen so that 
the /tm-sized grains in the extended halo (discussed 
next), assumed to originate from this cold compo- 
nent, are not too warm as to make the emission at 
25-35 /im inconsistent with the observations. The 



grain sizes range from 10 to 1000 /im. Large grains 
are assumed to exist because this is the location of 
the coUisional cascade that populates the extended 
halo discussed below. The estimated dust mass of 
this component is 0.12 Mq. This component is 
needed to fit the SED > 30 /im. It dominates the 
unresolved emission at 70 /xm and constitutes <9% 
of the total (unresolved-|-resolved) emission at 24 
/im. This component also accounts for the unre- 
solved sub-mm emission in Williams & Andrews 
(2006). 

• Extended halo. 
At 24 /im the cold component discussed above 
would be barely resolved, and a 70 /im it would 
be unresolved, therefore, it is inferred that most 
of the resolved 24 /tm and 70 /im emission comes 
from an extended halo component. It is modeled 
with a surface density, E ~ t~^ , characteristic of 
an unbound disk. The limited spatial resolution at 
70 /im does not allow determination of the inner 
radius; it is placed at Ri„ — 300 AU because it is 
assumed the grains in the halo originate in a coUi- 
sional cascade in the the unresolved cold disk. Rout 
is set to 1000 AU because the disk can be traced 
even beyond this distance (the best-fit outer radius 
ranges from 900 AU to 1800 AU). The grain sizes 
range from 1 to 10 /im (grains larger than the blow- 
out size are included to account for the effects of 
porosity on the response of the grains to radiation 
pressure). Assuming a grain size distribution fol- 
lowing n{b) ex &^''-^ and a maximum grain size of 
10 /xm, the ratio of the 70 and 24 /im fiuxes of the 
extended component is consistent with a minimum 
grain size of < 2 /im (similar to the blow-out size 
for HR 8799). The estimated dust mass of this 
component is 1.9-10^^ M®. 

In Section 5.4 we discuss the above results from Su et 
al. (2009)^ in the context of the dynamical simulations 
in Section 4.4; these simulations consider the five 
possible planetary configurations listed in Table 2 and 
described in Section 2.4. 



4. DYNAMICAL MODELING 

4.1. HD 128311 

Using numerical simulations we have studied the ef- 
fect of short and long-term planetary perturbations on 
the orbital stability of the dust-producing planetesimals. 

* Rcidcmcister et aL (2009) also carried out an analysis of the 
SED of HR8799. They assume a flat surface density distribution 
(S ~ r"), astronomical silicate composition, and particle size dis- 
tribution of n{b) oc b~^'^ with bmin = 5/im and hmax = 1000/im. 
They inferred the presence of a warm component with Rout = 10— 
15 AU, and a cold component with Ri„ = 75-120 AU and Rout 
= 125—170 AU. We will take this result with caution because: (1) 
this fit docs not include the deep IRS and spatially resolved MIPS 
observations (published after Reidemeister et al. 2009 paper); (2) 
they noted that the calibration of the (shallower) IRS data used 
is uncertain; and (3) we noted that their color-corrections for the 
IRAS observations were done for a characteristic temperature of 
5000 K, instead of the much smaller dust temperature (that would 
lead to smaller correction factors). 
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Fig. 7. — Allowed parameter space for the planetesimals' orbital 
elements. Shaded area: regions where test particle orbits are un- 
stable due to planet-crossing (red) or overlapping first order mean 
motion resonances (grey). The large black dots represent the two 
planets. The connected crosses are the maximum eccentricity at- 
tained by test particles on initially circular orbits, resulting from 
a numerical integration lasting 1 Myr. A zoom in to the region 
interior to the planets is shown in the right panel. The zone inside 
~0.3 AU and the zone beyond ~ 4 AU appears stable for potential 
locations of planetesimal belts. 

Notice that with the orbital elements adopted in this pa- 
per for HD 128311 (Models Al and A2 in Table 2), the 
planets are just barely outside their own planet-crossing 
zones. Because of the proximity of the planets' orbits to 
each other and to the 2:1 resonance, the analytical the- 
ory for the secular perturbations cannot be applied (as 
it was done for HD 38529 in Moro-Martm et al. 2007b), 
and therefore we study the stability of the orbits using 
numerical integrations. Figure 7 corresponds to Fit Al 
and shows the maximum eccentricity attained by test 
particles on initially circular orbits, resulting from gravi- 
tational perturbations by the two planets for an integra- 
tion time of 1 Myr. We can exclude the presence of plan- 
etesimals in the regions shaded with black dots because 
those orbits are chaotic due to overlapping first-order 
mean motion resonances (Wisdom 1980); similarly, the 
regions shaded with red dots are not expected to have sta- 
ble planetesimals because they are planet-crossing eccen- 
tric orbits. The stability of the potential dust-producing 
planetesimal belt was also studied for the orbital solu- 
tion in Fit A2 (Figure 8). For this simulation, we used 
500 particles uniformly spaced between 0.5 AU and 100 
AU, on initially circular orbits co-planar with the plan- 
ets, with angular elements chosen randomly between 
and 27r. Particles were removed if they approached the 
star closer than 0.5 AU, or approached a planet closer 
than its Hill radius. The orbits were integrated for ^ 
100 Myr using the multiple time step symplectic method 
skeel-SyMBA (Duncan, Levison & Lee 1998). The evolu- 
tion of the planets' semi-major axes and eccentricities is 
shown in Figure 18. The potential locations of planetesi- 
mal belts are the regions where the orbits are stable and 
where the maximum eccentricity is low (^0.3), ensuring 
long planetesimals' lifetimes. Figures 7 and 8 show that 
there are two such regions: (a) interior to ^ 0.3 AU and 
(b) exterior to ^ 4 AU (for a > 4 AU the maximum 
eccentricity always remains below 0.05). 

4.2. HD 202206 

The stability of the potential dust-producing planetes- 
imals was studied using numerical simulations of test 
particles with the planetary parameters listed in Table 
2 (Fit Bl). The simulations were done for 500 particles 
uniformly spaced between 0.5 AU and 100 AU, with 
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Fig. 8. — Results from the dynamical simulation of 500 test 
particles in the HD 128311 planetary system (Fit A2). A zoom 
in to the inner system is shown in the left panels. The numerical 
integration lasted 100 Myr. Top: test particle's lifetimes. Bottom: 
allowed parameter space for the planetesimals' orbital elements, 
where the shaded areas indicate regions where test particle's orbits 
are unstable due to planet-crossing {striped area) or overlapping 
first order mean motion resonances {dotted area); the squared sym- 
bols show the maximum eccentricity attained by test particles on 
initially circular orbits. 

an integration time of 100 Myr. The evolution of the 
planets' semi-major axes and eccentricities during that 
time span is shown in Figure 18. The results shown in 
Figure 9 indicate that the test particle's orbits are stable 
beyond ~ 6 AU, and that their maximum eccentricity 
remains < 0.3 for all semi-major axis; this indicates that 
the dust-producing planetesimals are located beyond 6 
AU. 



4.3. HD 82943 

We have studied numerically the stability of the poten- 
tial dust-producing planetesimals in the three planetary 
configurations discussed in Section 2.3 (listed in Table 2). 
The simulations were done for 500 particles uniformly 
spaced between 0.5 AU and 100 AU, with an integration 
time of > 70 Myr. The evolution of the planets' semi- 
major axes and eccentricities is shown in Figure 18. The 
results in Figures 10-12 indicate that, in the three plan- 
etary configurations considered, the test particle orbits 
are stable beyond '^ 3 AU, with maximum eccentricities 
always < 0.1. Long-lived, dust-producing planetesimals 
could therefore be located anywhere beyond 3 AU. 

4.4. HR 8799 

We have studied the stability of the potential dust- 
producing planetesimals using numerical simulations of 
test particles in the five planetary configurations dis- 
cussed in Section 2.4 and listed in Table 2. The sim- 
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Fig. 9. — Same as Figure 8 for 500 test particles in the HD 
202206 planetary system (Fit Bl); the numerical integration lasted 
100 Myr. 



Fig. 11.— Same as Figure 8 but for Fit C2 of HD 82943; the 
numerical integration lasted ~ 72 Myr. 
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Fig. 10. — Same as Figure 8 for 500 test particles in the HD 
82943 planetary system (Fit CI); the numerical integration lasted 
~ 82 Myr. 

ulations consisted of 1500 test particles: 500 particles 
uniformly spaced between 2 AU and the semi-major axis 
of the inner- most planet, 500 particles located between 
the innermost and outermost planet, and 500 particles 
uniformly spaced between the semi-major axis of the out- 
ermost planet and 300 AU. This latter value is taken from 
the results in the SED analysis in Su et al. (2009) sum- 
marized in Section 3.4. We assumed that the planets and 
the dust-producing planetesimals formed out of a thin 
disk and are co-planar. Particles were removed if they 
approached the star closer than 2 AU, or approached a 



Fig. 12.— Same as Figure 8 but for Fit C3 of HD 82943); the 
numerical integration lasted ~ 68 Myr. 

planet closer than its Hill radius. The orbits were inte- 
grated for 160 Myr (the upper limit to the stellar age 
estimated by Marois et al. 2008). The semi- major axes 
and eccentricities of the planet orbits do not evolve sig- 
nificantly during that time span (see Figure 18). The 
results regarding the test particles are shown in Figures 
13-17. For all the planetary configurations considered, 
the region between the orbits of the planets is dynami- 
cally unstable. The test particles' orbits are stable in two 
regions: (a) < 12 AU for fits Dl and D2, and < 10 AU 
for fits D3, D4, D5; and (b) > 110 AU for fits Dl, D2, 
D5, and > 150 AU for fits D3, D4. Regarding the maxi- 
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Fig. 13. — Same as Figure 8 for 1500 test particles in the HR 
8799 planetary system (Fit Dl); the numerical integration lasted 
160 Myr. The very slight rising trend of the maximum eccentricity 
at large semi-major axis seen in Figures 13, 14, 15, 16 and 17 is 
a numerical artifact that arises because the test particle simula- 
tions for HR8799 were carried out by skeel-SyMBA which gives or- 
bital elements in stellar-centric osculating elements; this produces 
a spurious eccentricity in nearly circular orbits because the star 
is wobbling, and its wobble velocity becomes a larger and larger 
fraction of the orbital velocity of test particles at larger and larger 
distances. (This effect has been corrected in Figures 8-12, where 
the osculating keplerian elements of the test particles are calculate 
in the barycentric frame.) 
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Fig. 14.— Same as Figure 13 but for Fit D2 of HR 8799. 
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Fig. 15.— Same as Figure 13 but for Fit D3 of HR 8799. 
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Fig. 16.— Same as Figure 13 but for Fit D4 of HR 8799. 

5. DISCUSSION: POTENTIAL LOCATION OF THE 
DUST-PRODUCING PLANETESIMAL BELTS 

We now discuss the potential location of the dust- 
producing planetesimal belts in each for the four 
multiple-planet systems studied in this paper. 



mum eccentricity of the test particles in the dynamically 
stable regions, Figures 13-17 show that Cmax ^ 0.3, even 
though for some fits there are peaks at ~ 5 AU and ^170 
AU (due to secular resonances) . Given the low maximum 
eccentricities, these dynamically stable regions could be 
possible locations of long-lived, dust-producing planetes- 
imals (except at ~ 5 AU in fit D3 where the Cmax reaches 
0.4 and where the planetesimals lifetimes might be short 
under coUisional evolution). 



5.1. HD 128311 

We favor models that assume a single grain radius of 
10 /xm because they represent large blackbody grains in 
general, and there is no evidence of a significant popu- 
lation of small grains. With these models, the observed 
SED can be fitted by a dust disk with Rout — 10 AU 
and Ri„ > 5 AU, Rout = 50 AU and R„ > 2 AU, or 
Rout = 100 AU and Ri„ > 2 AU. The dynamical simula- 
tions help constrain further the possible location of the 
dust-producing planetesimals because the effects of the 
planets extend into the regions allowed by the SED mod- 
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Fig. 17.— Same as Figure 13 but for Fit D5 of HR 8799. 

els: they predict two stable niches where planetesimals 
could be long-lived, beyond ~ 4 AU and inside ^ 0.3 
AU. Regarding the latter, to account for the lack of ex- 
cess emission at A < 33 ^m, the SED modeling excludes 
the presence of a significant population of dust-producing 
planetesimals in this region. To set tighter constraints to 
the location of the planetesimals in HD 128311 there is 
the need to obtain spatially resolved images and/or accu- 
rate photometric points in the 33 /zm-70 /zm range and 
in the sub- mm. As seen in Figures 2 and 3, observa- 
tions with the recently launched Herschel/PACS would 
be very valuable for this purpose. 

5.2. HD 202206 

The dynamical model indicates that planetesimals or- 
bits are stable beyond ^ 6 AU, and that their maximum 
eccentricity remains ;$ 0.1 for all semi-major axes, i.e. 
that planetesimals could be long-lived beyond 6 AU. This 
result agrees with the conclusions from the SED model- 
ing. These models, which are well-constrained because 
of the presence of a small excess beyond 25 /im, result in 
a relatively narrow range of planetesimal belts that can 
fit the observations, ranging from a 50 AU disk with an 
inner cavity 10-20 AU in size, to a 100 AU disk with an 
inner cavity 6-10 AU in size. We conclude that the grav- 
itational perturbations of the detected planets might be 
responsible for the inner edge of the dust disk. 

5.3. HD 82943 

The observed SED, lacking emission at A < 33 /im, can 
be fitted by a dust disk composed of single grains 10 fiia 
in size, with an inner cavity with Ri„ >16 AU for a disk 
with Rout = 50 AU, or R„i = 12 -26 AU for a disk with 
R-out = 100 AU (a compact 10 AU dust disk is excluded). 
Because the dynamical modeling of test particles in the 
three planetary configurations considered suggests that 
planetesimals could be stable and long-lived beyond '^ 
3 AU, we conclude that the gravitational perturbations 
from the planets (located within 2.1 AU of the star) do 
not extend far enough to have a significant effect to desta- 



bilize any debris system that is seen in the infrared excess 
emission. 

5.4. HR 8799 

For the five planetary configurations considered, the 
dynamical modeling of test particles suggests that plan- 
etesimals could be stable and long-lived at semi-major 
axes < 10-12AU and > 110-150 AU. This helps to fur- 
ther constrain the dust disk solutions based on the SED 
fitting by Su et al. (2009), who proposed: (1) the pres- 
ence of an unresolved warm disk 15 AU in size, if the 
grain sizes range from 1 to 4.5 ^m, or a more compact 
10 AU disk, if all the grains are bound (> 2 /im); and 
(2) the presence of an unresolved cold disk with an inner 
edge at ^ 90 AU. The dynamical models favor an inner 
dust disk component of bound grains and Kout ^ 10 AU 
and a colder component with Ri„ ~ 110-150 AU (rather 
than ~ 90 AU). Regarding the cold component, we find 
there is no tension between Ri„ derived from the dynam- 
ical models {^ 110-150 AU) and that derived from the 
analysis of the SED and the surface brightness radial pro- 
files. Even though Su et al. (2009) favors a value of Ri„ 
~ 90 AU, Figure 19 shows that the observations can also 
be fitted with Ri„ = 110 AU (increasing the total dust 
mass) and Ri„ ~ 150 AU (changing the minimum grain 
size from 10 /im to 8 /im); this is not surprising because 
all these inner radii are smaller than the MIPS 70 /im 
pixel size. 

The presence of the outer planetesimal disk may help 
constrain planet formation scenarios that have been pro- 
posed for the HR8799 system. One scenario is planet- 
planet scattering: in addition to the difficulty to result 
in a stable system with low eccentric planets (Dodson- 
Robinson et al. 2009), this model may not be able to ac- 
count for the presence of the outer planetesimal disk. A 
second scenario is long-ranged outward migration in reso- 
nance (Crida et al. 2009): future high-resolution imaging 
of the cold component of the HR 8799 debris disk may be 
able to provide evidence of dust-producing planetesimals 
trapped in MMRs with the outermost planet, a signpost 
of outward planet migration. 

Table 3 and Figure 1 summarize the possible planet- 
planctesimals configurations of HD 128311, HD 202206, 
HD 82943 and HR 8799, compared to that of the other 
three multi-planet systems known to harbor dust - HD 
38529, HD 69830 and the Sun. In some cases, the SED al- 
lows for the presence of both narrow and wide belts. We 
favor the latter because the narrowest belts might only 
be justified in the presence of additional planetary per- 
turbers for which we have no evidence so far. However, 
one should keep in mind the degeneracy: to set tighter 
constraints to the location of the planetesimals there is 
the need to obtain spatially resolved images and/or accu- 
rate photometric points in the 30 /xm-70 fim range and 
in the sub- mm. Observations with Herschel/PACS and 
ALMA will be very valuable for this purpose. 

6. CONCLUSIONS 

In this paper we have studied the possible planet- 
planctcsimal configurations of four multi-planet systems, 
of which three are radial-velocity systems - HD 128311, 
HD 202206, HD 82943 - and one' is a directly imaged 
system - HR 8799. We have quantified where the zone 
of influence lies of planets on the dust-producing plan- 
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Fig. 18. — Long-term evolution of the planets' semi-major axes and eccentricities for the planetary configurations in Table 2. From 
top to bottom, left to right the panels correspond to fits A2, Bl, C1-C3, D1-D5, for HD 128311, HD 202206, HD 82943 and HR 8799, 
respectively. The models were run for 100 Myr (HD 128311), 100 Myr (HD 202206), >70 Myr (HD 82943) and 160 Myr (HR 8799). The 
colors correspond to planet b (black), c (orange) and d (blue). 



etesimals: for HR 8799 it extends to nearly 20 AU from 
the orbit of the outermost planet, while for the three 
radial- velocity systems it extends to about 4 AU. A previ- 
ous paper that studied HD38529, another radial- velocity 
multi-planet system, found that the influence of the plan- 
ets in this case extends out to ~ 10 AU (determining the 
inner edge of the disk) , and becomes dominant again at 
a ~ 55 AU (due to a secular resonance that probably 
determines the outer edge of the dust disk (Moro-Martin 
et al. 2007b). We conclude that radial- velocity multi- 
planet systems generally have zones of influence within a 
few to ten AU; more precise determinations will require 
individual modeling of a system. The influence can be 
extended much further through secular resonances and 
similar behavior. 

For the three radial- velocity multi-planet systems stud- 
ied in this paper, we have constructed fits to the spectral 
energy distributions of the debris disks. If we use as- 
tronomical silicates with a size distribution down to the 
blow-out size, the emitting regions of the disks must be 
so far from the star that they are well outside the zones 
of influence of the radial velocity planets. This behavior 
depends critically on the optical properties of the grains 
(size distribution and optical constants of the grain ma- 
terial). If we were to adopt optical constants typical of an 
ice-silicate mixture instead of astronomical silicates, the 



dust may be located closer to the star. In all three cases, 
we find that single-size 10 /im astronomical silicates re- 
produce well the observed SED, and could lie at the edge 
of the zone of influence of the planets and therefore the 
disk would be sculpted by them. 
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AU, Rout = 1000 AU, a surface density given by S ~ r ^, M^^^j = 1.! 
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TABLE 3 

Stars with evidence of harboring a multiplanet system and 
planetesimals 



Source 



SpTyp 
Age (Gyr) 



Excess A (fim) 
Ld„,t/L.(10-*) 



Planetesimals' Location^ (A-U) 
lOfira grains size distribution 





Planetary Configuration'' 


Fit 1 


Fit 2 


Planet 


Planet 


c d 


e b c d 



Fit 3 
Planet 
c d 



Sun'' 



HD 128311 KO V 70 >5 10== >52 100== M 2.19 3.22 

0.5 0.13-0.27 >2 50-= a 1.10 1.76 ■•■ 

>2 100° e 0.25 0.17 


1.61 3.18 
■ 1.11 1.73 
• 0.36 0.21 


HD 202206 G6 V 25-70 10-20 50 M 16.6 2.19 ■■• 
5.6 6-10 100 a 0.80 2.51 ■•■ 

e 0.44 0.07 




HD 82943 GO V 70 16 50° M 1.46 1.73 ■ ■ • 
5 0.88-1 12-26 100 a 0.75 1.19 

e 0.45 0.27 ■•■ 


■ 1.68 1.87 0.49 1.70 1.74 0.35 
0.75 1.20 2.12 0.74 1.20 1.19 

■ 0.39 0.11 0.018 0.36 0.19 0.078 


HR 8799° A5 V 8-850 warm: 6 to 10 M 7 10 10 
0.03-0.16 2.3 cold: 110-150 to 300 a 67.1 38.0 23.4 

e 0.09 ■ 


8.02 11.9 8.89 9.71 7.96 7.40 

• 68.4 39.6 24.2 67.7 31.0 30.8 

0.01 0.01 0.07 0.01 0.25 0.27 


HD 38529 G8 III/IV 70 15 50 M 0.85 13.2 ■ ■ • 
3.5 0.36 8 100 a 0.13 3.74 ■•■ 

Rsub 500 e 0.25 0.35 ■•■ 




HD 69830 KOV 8-35 0.93-1.16 M 0.03 0.04 0.06 ■ 
4-10 2 a 0.08 0.19 0.63 

e 0.10 0.13 0.07 • 





G2 V 

4.57 



(0.1-l)-10-3 
(l-10)-10-3 



2-4 (Asteroid Belt) 
35-50 (Kuiper Belt) 



M 



1 

5.20 

0.05 



0.30 0.05 0.05 
9.58 19.2 30.1 
0.05 0.04 0.01 



"" Inferred planetesimals' location. HD 128311: this work; HD 202206: this work; HD 82943: this work; HR 8799 : Su ct al. (2009) and this work; HD 38529: Moro-Martin et al. 
2007b; HD 69830: Lisse et al. 2007. 

^ Planetary parameters: M, a and e are planet mass in Mj„p, semi-major axis in AU and eccentricity, respectively. The rest of the orbital elements are listed in Table 2. 

° For HR 8799, this table only includes the planetary configurations where the planets have their nominal masses; the SED reveals the presence of a warm component and a cold 
component of the dust disk (from Su et al. 2009). 

'^ Fractional luminosities of asteroidal dust and Kuiper belt dust from Dermott et al. (2002) and Stern (1996). 

° Narrow belts of width 10% of the disk radius give valid fits to the SED (but see discussion in Section 3.1). 



